clear;clc;close all;
%% data (from PHYSICAL REVIEW A 67, 022505 (2003))
load xi_Delta_data.mat

s0=1.00624;
lambda0=exp(pi/s0);
ainv_disc=exp(Delta_disc/2/s0).*cos(xi_disc);

hbar=1;
m=1;
xi_star=1;

Energy_disc=@(n) -ainv_disc.^2+xi_star^2*lambda0^(-2*n)*exp(Delta_disc/s0);

%% fitting results (the fitting process is done by Curve Fitting Toolbox 23.2)
load fitting_data.mat

xi1=-pi:.005:-5*pi/8;
xi2=-5*pi/8:.005:-3*pi/8;
xi3=-3*pi/8:.005:-pi/4;
xi=[xi1 xi2 xi3];

z=(xi1+pi).^2.*exp(-1.*(xi1+pi).^(-2));
y=xi2+pi/2;
x=sqrt(-xi3-pi/4);

Delta1=a2*z.^2+a1*z+a0;
Delta2=b0+b1*y+b2*y.^2+b3*y.^3;
Delta3=c0+c1*x+c2*x.^2;
Delta=[Delta1 Delta2 Delta3];

ainv=exp(Delta/2/s0).*cos(xi);
Energy=@(n) -ainv.^2+xi_star^2*lambda0^(-2*n)*exp(Delta/s0);

%% plotting 1
figure
hold on;box on;grid on
plot([0 0.5],[0 -0.5],"-.",LineWidth=2.5,Color="#808080")
plot(sign(ainv).*sqrt(abs(ainv)),-(Energy(0)).^0.25,'-',LineWidth=3.5,Color="#4169E1")
plot(sign(ainv_disc).*sqrt(abs(ainv_disc)),-(Energy_disc(0)).^.25,"o",MarkerSize=8,LineWidth=2.5,Color="#00008B")
xlabel('sgn$(a)(\kappa_*|a|)^{-1/2}$',Interpreter='latex');
ylabel('$-(mB_3/\kappa_*)^{1/4}$',Interpreter='latex');
set(gca,'LineWidth',2.5);
set(gca,'FontSize',20);
xlim([-0.2 0.5])
lgd=legend(["","fit","data"]);
set(lgd,FontSize=20,location="southwest");
set(gcf,'position',[300,100,500,400])

%% plotting 2
s=2.2;
figure
hold on;box on
plot([-0.5 0.5],[0 0],"-",LineWidth=2.5,Color="#808080")
plot([0 0],[-0.5 0.5],"-",LineWidth=2.5,Color="#808080")
plot([0 0.5],[0 -0.5],"-.",LineWidth=2.5,Color="#808080")
plot(sign(ainv).*sqrt(abs(ainv)),-(Energy(0)).^0.25,'-',LineWidth=3.5,Color="#00008B")
plot(sign(ainv).*sqrt(abs(ainv))/s,-(Energy(0)).^0.25/s,'-',LineWidth=3.5,Color="#4169E1")
plot(sign(ainv).*sqrt(abs(ainv))/s^2,-(Energy(0)).^0.25/s^2,'-',LineWidth=3.5,Color="#87CEFA")
plot(sign(ainv).*sqrt(abs(ainv))/s^3,-(Energy(0)).^0.25/s^3,'-',LineWidth=3.5,Color="#ADD8E6")
plot(sign(ainv).*sqrt(abs(ainv))/s^4,-(Energy(0)).^0.25/s^4,'-',LineWidth=3.5,Color="#E6E6FA")
xlabel('sgn$(a)(\kappa_*|a|)^{-1/2}$',Interpreter='latex');
ylabel('$-(mB_3/\kappa_*)^{1/4}$',Interpreter='latex');
set(gca,'LineWidth',2.5);
set(gca,'FontSize',20);
xlim([-0.11 0.11])
ylim([-0.16 0.04])
yticks([-0.16 -0.12 -0.08 -0.04 0 0.04])
set(gcf,'position',[300,80,880,640])

%% plotting 3
s=2.2;
figure
hold on;box on
plot([-0.5 0.5],[0 0],"-",LineWidth=2.5,Color="#808080")
plot([-0.5 0.5],[0 0],"-",LineWidth=2.5,Color="#808080")
plot([0 0],[-0.5 0.5],"-",LineWidth=2.5,Color="#808080")
plot([0 0.5],[0 -0.5],"-.",LineWidth=2.5,Color="#808080")
plot([0 0.5],[0 -0.6],"-.",LineWidth=2.5,Color="#808080")
X1=sign(ainv).*sqrt(abs(ainv));
Y1=-(Energy(0)).^0.25;
plot(sign(ainv).*sqrt(abs(ainv)),-(Energy(0)).^0.25,'-',LineWidth=3,Color="#483D8B")
plot(X1(1:438)/0.94,Y1(1:438)/0.94,'-',LineWidth=3,Color="#DC143C")
plot(X1(1:438)/0.78,Y1(1:438)/0.78,'-',LineWidth=3,Color="#FF8C00")

X2=sign(ainv).*sqrt(abs(ainv))/s;
Y2=-(Energy(0)).^0.25/s;
plot(sign(ainv).*sqrt(abs(ainv))/s,-(Energy(0)).^0.25/s,'-',LineWidth=3,Color="#483D8B")
plot(X2(1:438)/0.94,Y2(1:438)/0.94,'-',LineWidth=3,Color="#DC143C")
plot(X2(1:438)/0.78,Y2(1:438)/0.78,'-',LineWidth=3,Color="#FF8C00")

X3=sign(ainv).*sqrt(abs(ainv))/s^2;
Y3=-(Energy(0)).^0.25/s^2;
plot(sign(ainv).*sqrt(abs(ainv))/s^2,-(Energy(0)).^0.25/s^2,'-',LineWidth=3,Color="#483D8B")
plot(X3(1:445)/0.94,Y3(1:445)/0.9,'-',LineWidth=3,Color="#DC143C")
plot(X3(1:456)/0.78,Y3(1:456)/0.7,'-',LineWidth=3,Color="#FF8C00")


xlabel('sgn$(a)(\kappa_*|a|)^{-1/2}$',Interpreter='latex');
ylabel('$-(mB_3/\kappa_*)^{1/4}$',Interpreter='latex');
set(gca,'LineWidth',2.5);
set(gca,'FontSize',20);
xlim([-0.135 0.135])
ylim([-0.18 0.04])
yticks([-0.16 -0.12 -0.08 -0.04 0 0.04])
set(gcf,'position',[300,80,880,640])
